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The properties of inhomogeneous nuclear matter are investigated considering the 

■ self-consistent Skyrme Hartree-Fock approach with inclusion of pairing correlations. 
For a comparison we also consider a relativistic mean field approach. The inhomo- 
geneous infinite matter is described in terms of cubic Wigner-Seitz cells, which leads 
to a smooth transition to the limit of homogeneous nuclear matter. The possible 

Q_i! existence of various structures in the so-called pasta phase is investigated within this 

| self-consistent approach and a comparison is made to results obtained within the 

' Thomas- Fermi approximation. Results for the proton abundances and the pairing 

properties are discussed for densities for which clustering phenomena are obtained. 
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I. INTRODUCTION 

The crust of neutron stars is a very intriguing object for theoretical nuclear structure 

■ physics, as it contains the transition from stable nuclei in the outer crust to a system of 
homogeneous nuclear matter, consisting of protons, neutrons and leptons in /9-equilibrium, 
in the inner part of this crust. The question of how matter consisting of isolated nuclei melts 
into uniform matter with increasing density has evoked a large number of studies[l|, 0, 0, 0]. 

■^J- ■ Already at moderate densities the Fermi energy of the electron is so high that the /9-stability 
enhances the neutron fraction of the baryons so much that a part of these neutrons drip 
out of the nuclei. This leads to a structure, in which quasi-nuclei, clusters of protons and 
neutrons, are embedded in a sea of neutrons. In order to minimise the Coulomb repulsion 
between the protons, the quasi-nuclei form a lattice. 

Therefore one typically describes these structures in form of the Wigner-Seitz (WS) cell 
approximation. One assumes a geometrical shape for the quasi-nuclei and determines the 
nuclear contribution to the energy of such a WS cell from a phenomenological energy-density 
functional. Such Thomas-Fermi calculations yield a variety of structures: Spherical quasi- 
nuclei, which are favoured at small densities, merge with increasing density to strings, which 
then may cluster to parallel plates and so on. These geometrical structures have been the 
origin for the popular name of this phase: Pasta phase. 

Such Thomas-Fermi calculations, however, are very sensitive to the surface tension under 
consideration. Furthermore they do not account for characteristic features of the structure 
of finite nuclei, like the shell-effects. Shell effects favour the formation of closed shell systems 
and may have a significant effect on the formation of inhomogeneous nuclear structures in 
the crust of neutron stars. These shell effects are incorporated in self-consistent Hartree-Fock 
or mean field calculations, which can treat finite nuclei, infinite matter and inhomogeneous 
structures in between within a consistent frame based on an effective nucieon-nucleon inter- 
action. Such calculations employing the density-dependent Skyrme forces [1, @] have been 
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done more than 25 years ago by Bonche and Vautherinj3] and by a few other groups. 

These studies show indeed that shell effects have a significant influence on details like the 
proton fraction of the baryonic matter in the inhomogeneous phase [8|. They also provide 
the basis for a microscopic investigation of properties beyond the equation of state. This 
includes the study of pairing phenomena, excitation modes and response functions as well 
as the effects of finite temperature. 

Self-consistent Hartree-Fock calculations for such inhomogeneous nuclear structures have 
typically been performed assuming a WS cell of spherical shape. This assumption of quasi- 
nuclei with spherical symmetry reduces the numerical work-load considerably. However, 
it does not allow the exploration of quasi-nuclear clusters in form of strings or plates as 
predicted from Thomas-Fermi calculations. Furthermore the limit of homogeneous matter 
can not be described in a satisfactory manner in such a spherical WS cell. Employing the 
representation of plane wave single-particle states in terms of spherical Bessel functions leads 
to a density profile, which, depending on the boundary condition chosen, exhibits either a 
minimum or a maximum at the boundary of the cell. Bounche and Vautherin0] therefore 
suggested to use a mixed basis, for which, depending on the angular momentum, different 
boundary conditions were considered. However, even this optimised choice leads to density 
profiles with fluctuations (§]. 

Therefore the investigations presented here consider cubic WS cells, which allows for the 
description of non-spherical quasi-nuclear structures and contains the limit of homogeneous 
matter in a natural way. Self-consistent Hartree-Fock calculations are performed for /3-stable 
matter at densities for which the quasi-nuclear structures discussed above are expected. For 
the nuclear Hamiltonian we consider various Skyrme forces but also perform calculations 
within the effective relativistic mean-field approximation. Special attention will be paid to 
the comparison between results obtained in the Hartree-Fock approach and corresponding 
Thomas-Fermi calculations. 

After this introduction we will briefly review the Hartree-Fock approximation using 
Skyrme interactions and the technique used to solve the equations resulting from this ap- 
proach employing the imaginary time step method in section 2. We then turn to the rela- 
tivistic mean field approach and the adaption of the imaginary time step method to be used 
within this relativistic framework. After a short description on the inclusion of pairing cor- 
relations in section 4, we present results in section 5. The main conclusions are summarised 
in the final section 6. 



II. SKYRME HARTREE FOCK CALCULATIONS 



A. Energy Functional 

The Skyrme-Hartree-Fock approach has frequently been described in the literature 0, @, 
0, H|. Therefore we will restrict the presentation here to a few basic equations, which will 
define the nomenclature. The Skyrme model is defined in terms of an energy density TC(r), 
which can be split into various contributions 0, [l(| 

T~t = TLk + T~io + 7^3 + Tt c tt + 7ifi n + H so + Hcoui, (1) 

where TCk is the kinetic energy term, Ho a zero range term, 7i 3 a density dependent term, 
7i c s an effective mass term, Hfi n a finite range term and 7i so a spin-orbit term. These terms 
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are given by 



2m 

Ho = \t [(2 + x )p 2 -(2x + l)(p 2 p + p 2 n )], 
H 3 = ±t 3 p a [(2 + x 3 )p 2 -(2x 3 + l)(p 2 p + pl 
n cS = l[t 1 (2 + x l )+t 2 (2 + x 2 )]rp 

+ | [h(2x2 + 1) - ti(2xi + 1)] [r p p p + r n p n ] , 
ftfin = ~^[3t 1 (2 + x 1 )-t 2 (2 + x 2 )]pAp 

+^ [3t 1 (2x 1 + 1) + t 2 {2x 2 + 1)] [p p Ap p + p n Ap n ] , 

n so = -\W Q [ P VJ + Pp VJ p + p n VJ n ]. (2) 



The coefficients ti, X{, Wo, and a are the parameters of a generalised Skyrme force [11]. The 
energy density of eq.([T]) contains furthermore the contribution of the Coulomb force, Hcoui, 
which is calculated from the charge density pc as 

^coui = yPc(r) y d r |r _ r/| - — [-) Pc - ( 3 ) 

Here the exchange part of the Coulomb term is calculated within the Slater approximation. 
Following [lH the center-of-mass recoil energy has been approximated as — Y^p 2 /2Am. 

The densities p, r, and J are defined in terms of the corresponding densities for protons 
and neutrons p = p p +p n , i~ p +T n and J = J p + J n . If we identify the isospin label (q = n,p), 
the corresponding matter densities are given by 

P,W = E V q M(r,s)\ 2 , (4) 

k,s 

where (p g k (r,s) is the single-particle wave function with orbital, spin and isospin quantum 
numbers k, s and q. The occupation factors rfj. are determined by the Fermi energy and 
the desired scheme of occupation (see discussion below). The kinetic energy and spin-orbit 
densities are defined by 



r 9 (r) = 4\V<pl(r,s)\ 2 , (5) 

k,s 

J q (r) = -i £ 4 (<pl)*{r, s') V^(r, s) x ( s '|<r | s ) . (6) 



The gradient of the spin-orbit density VJ = VJ P + VJ n can be directly evaluated without 
first calculating J: 

VJ 9 (r) = -i J2 Vl VW(r, s') x V^(r, s) ■ (s>|s>. (7) 



We left out the spin-gradient term [lOj, which is cumbersome to evaluate in three- 
dimensional calculations numerically and not very important. 



4 



The single-particle wave functions are determined as solutions of the Hartree-Fock equa- 
tions 

{~ v 2^y v + Uq{r) ~ 1 Wq{r) ' (v x a) ] ^ (r) = el ^ (r> s) ' (8) 

with an effective mass term m*(r), which depends on the H c s part of the energy density 
functional 

h 2 ft 2 

— + ±[t 1 (2 + x 1 )+t 2 (2 + x 2 )]p(r) 



2m* (r) 2m 



+ |[t 2 (l + 2x 2 ) -h(l + 2x x )\ p q (r), (9) 



a nuclear central Potential 



U q (r) = |t [(2 + x )p - (1 + 2x )p g ] 
+ ^t 3 (2 + x 3 )(2 + a)p Q + 1 
- ±t 3 (2x 3 + 1) [2p> + ap"- 1 (pi + pi)] 
+ \[t 1 {2 + x 1 ) +t 2 {2 + x 2 )] r 
+ l [i 2 (2x 2 + l)-t 1 (2x 1 + l)] r q 
+ ±[h(2 + xi)-3t 1 (2 + x 1 )] Ap 
+ ^ [3ti(2xi + 1) + t2(2a; 2 + 1)] Ap q 
-\W [VJ + VJ q ] 

+ <J 9tP Vboui (10) 



with the Coulomb field 



- e 2 f-V /3 p^ 3 (ID 



7r y 

and a spin-orbit field: 

W q (r)= l Wo (Vp + V Pq ) (12) 

B. Imaginary Time Step 

Various different methods have been developed to solve the Hartree-Fock equations. Fre- 
quently the single-particle wave functions are expanded in a basis like e.g. the eigenfunctions 
of an appropriate harmonic oscillator. This is appropriate for describing the wave functions 
for single-particle states, which are deeply bound. It is not so appropriate for the description 
of weakly bound or unbound single-particle states, since the asymptotic behaviour of the 
harmonic oscillator basis states is not appropriate for these states. 

This can be cured by employing the eigenstates of a spherical box with an appropriate 
radius which can also be considered as a Wigner Seitz cell for describing periodic 

systems. Such a spherical box, however, is not appropriate for the description of deformed 
nuclei and nuclear structures as they are expected for the pasta phase in the crust of neutron 
stars. This, as well as the problems with the boundary conditions in a spherical WS cell 
discussed already in the introduction calls for a cartesian WS cell. 
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The single-particle wave functions in such a cartesian WS cell can be represented by its 
values on a discretized mesh in this cell. The spacings between the mesh points, Ax, Ay 
and Az correspond to truncations in momentum space. Smaller values for these spacings 
account for larger momentum components in the wave functions. The obvious disadvantage 
of such calculations is the huge amount of mesh points which has to be taken into account. 
Therefore one needs a fast iterative procedure for the solution of the self-consistent Hartree- 
Fock equations, which evaluates only the desired states. 

Davies et al. presented in [12|] an efficient method for this problem, the imaginary time 
step method, which we want to outline briefly. The origin for the name of this method is the 
analogy to the time-dependent Hartree-Fock (TDHF) method which solves the equations 

ih^- = H(t)(p k (t), k = l,...,A (13) 

for an orthonormal set of A wave functions {</?&}, arid a Hamiltonian H which depends on the 
time t. This is the case when we identify H(t) e.g. with the Hartree-Fock (HF) Hamiltonian 
represented in eq.flBJ), which depends on t as it depends on the resulting wave-function (pk(t) 
in a self-consistent way. These equations are discretized in time introducing a time step 
At, with t n = nAt. Then the time evolution of the set of wave functions {fk} ma y be 
approximated by the iterative procedure 

k n+1) > = exp (-1 At \<p£>), k = 1, . . . , A (14) 

in which yj^ represents the wave-function (p^ at the time t n and if ( n+ s) denotes the numer- 
ical approximation to the Hamiltonian Hit) at the time (n + \)At. The idea of Davies et al. 
was to replace the time step At by the imaginary quantity —iAt. Introducing the positive 
parameter A = At/H the procedure for the imaginary time step gets 

|y3i n+1) > = exp (-A |^ n) >, k=l,...,A (15) 

where {<Pk L+1 ^} is not any more an orthonormal set of wave functions since the imaginary time 
operator exp ( — XH^ n+ 2n is not unitary. Applying the Gram-Schmidt orthonormalization 
method O we get the orthonormal set by 

= 0\(pt +l) ) k = l,...,A. (16) 

This procedure converges leading to those eigenfunctions of the Hamiltonian H, which cor- 
respond to the lowest A eigenvalues of the Hamiltonian H . 

In practical applications the Hamiltonian H^ n+ 2> is replaced by the Hamiltonian H^ n ' of 
the n-th step, which makes the calculation fast keeping the algorithm stable. After this 
replacement, we get the following operation on the wave functions 

^k +1) = O ( exp ( - A ffW) tp™ ) k — l,...,A, (17) 

For numerical application one has to truncate the exponential series to a certain order. In 
earlier HF calculations the gradient method was used with the operation 0(1 — XH) on the 
wave functions ((sj). If one truncates the exponential series in the imaginary time step going 
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beyond the first order one obtains an improvement of the gradient method. Davies et al. 
recommended a truncation to 4th or 5th order for a HF calculation of 40 Ca together with a 
time step At = 4.0 x 10 -24 s and a mesh size of 1.0 fm. 

In our calculations we used the same mesh size as Davies et al. but the convergence 
got worse since we consider in our studies a larger number of nucleons, which implies a 
larger number of wave functions A have to be evolved. Hence we truncated the exponential 
operator at 9th order and the time step At was set to 2.0 x 10~ 24 s. For the check of 
convergence the mean square deviation of the single particle energies for N Nucleons and r\ k 
the occupation probability is calculated by 

AlfM = (jf: m ((^\H^W) - 2 (18) 

which provides a better criterion as calculating energy differences. 

The HF equations have been solved by discretization in coordinate space within a cubic 



Wigner-Seitz cell similar to pjj with . periodic boundary conditions. The box sizes typically 
considered vary from 2 x 10 fm to 2 x 16 fm. This technique is able to allow for general 
deformations of the quasi-nuclear structures. The densities we are considering requires to 
account for around 1500 nucleons, which implies that up to A = 2600 wave functions had 
to be evolved to account for pairing correlations with occupation probabilities t] k different 
from zero. 



To decrease the numerical effort we assume two symmetries like in [11] : 

• parity 

PLp k (r,s) = ip k (-r,s) = Pk Lp k (r, s), p k = ±1; (19) 

• z-signature 

exp{in(J z - \)}ip k (x, y, z, s) = <r<p k (-x, -y, z, s) ^ 

= T) k ip k (x,y,z,s), T] k = ±l. 

These symmetries still allow triaxial deformations and reduce the calculation to the pos- 
itive coordinates in each direction. As additional symmetry time-reversal-invariance is 
assumed for the time-reversed pairs ip k , and ip k : 

<p- k (r, s) = (f<p k )(r, s) = a<pl(r, -s). (21) 

Summarising this symmetries it is sufficient to solve the HF equations for one wave function 
of the time-reversed pairs. We choose the positive z-signature orbital for which we get the 
symmetries summarised in table [B The wave functions ^(r, s) are realized as complex Pauli 
spinors. The reflections at the x = and y = planes are realized by the parity operator of 
the real part together with complex conjugation. 

The iteration is performed with accurate numerical methods. For the differential opera- 
tors 11-point formulas are used, which have been derived by eliminating errors for functions 
/ with f(x) = x n up to a certain no G N. The ansatz for the numerical approximation of 
the derivatives on an equidistant mesh with the points and fa = f(xi) is for the first 
derivative 

s'<*> M (s) /w = t^2]k («) 

\ / num j = i J 
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x=0 y=0 z=0 


Re ip k (r,+l) 
Im (p k (r, +i) 
Re y k {r,-\) 
Im <^ fc (r, -i) 


+ + Pfe 

- - Pfc 

- + ~Pk 
+ ~ ~Pk 



TABLE I: Parity properties of the Pauli spinors with respect to the coordinate planes 

with N = 5 for 11-point formula and aj the coefficients of the formula. Requiring that the 
approximation gets equal up to a certain no G N we obtain a linear equation. Inserting the 
result in the ansatz we finally obtain 



d_ 

dx 



f(xi) (23) 



= ^(tM6o(11/*+5 - 4500/ i+2 + 16350/ i+1 

- 16350/ i _ 1 + 4500/^2 -ll/i_ 5 ) 
+ y^(-445/, +4 + 2950/.+3 - 2950/^3 + 445/U))- 

For the second derivative used in the laplacian the ansatz is 

N 



d_ 

dx 



(24) 



and finally the formula gets 
d 2 



dx 2 



f{xi) (25) 



I1UII1 



(Ax) 



L-^(ll/, +5 - 11250/, +2 + 81750/ m 



+ 81750/^ - 11250/U + ll/i_ 5 ) - S/' 
+ 33^4s(- 1335 ^+4 + 11800/ i+3 + 118002950/ J _3 - 1335/,_ 4 )). 



In the case of the Wigner-Seitz cell calculations charge neutrality is assumed and electrons 
are taken into account as relativistic Fermi gas which contribute to the charge density 
Pc(r) = p p (r) — p e . For the calculation of finite nuclei the electrons are not taken into 
account (p e = 0). 

There are different methods to solve the Poisson equation 

- AV c (r) =47r Pc (r). (26) 

It turned out that the numerically most accurate and stable method is the integration 
applying the Green's function for this problem 

Vc(r)= I dr'" Pc {r') —^—r (27) 
Jv \r — r \ 
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Unfortunately, this integral has lots of singularities, but it can be rewritten. First, the 
Green's function is written as (fl3|): 

' =lA r ,\r-r'\. (28) 



\r — r 



Then the integral is transformed by Green's theorem for scalar functions / and g defined on 
a Volume V with closed surface A = dV (UA 



[ dV(fAg)- [ dV(gAf)= I dA-(fVg-gVf). (29) 

JV JV JA=dV 

Identifying and / = pc( r ') and g = -Ar — r'\ the final result gets 

V c {r) = \ [ dr' 3 Ap c (r')\r-r'\ 

2 J \ , (30) 
+ - f dA (pc(r') V r '\r - r'\ - \r - r'\ V r >p c (r')) , 

* JA=dV 

which has no singularities. Altogether the result of this transformation behaves very well 
in numerical calculations and the numerical result is practically the same as the exact one. 
For finite nuclei it is possible to drop the boundary integrals as has already been discussed 
by Vautherin fl~3l ]. 

We tested the computer program for the parameter set Skyrme III by comparing results 
for finite nuclei with those of 11] . Additional tests have been performed using the parameter 
set SLy4 0- 



III. RELATIVISTIC MEAN FIELD CALCULATIONS 

In order to test the sensitivity of the results on the model under consideration we also in- 
vestigated the quasi-nuclear structures in the crust of neutron stars employing the relativistic 
mean field approach in a cubic box. 



A. From the Lagrangian to the Dirac Equation 



The relativistic mean field approach is based on a Lagrangian is similar to that in [15 
and consists of three parts: Lagrangian for the free baryons Cb, the free mesons Cm and 
the interaction Lagrangian £j nt : 

C = Cb + Cm + C int , (31) 

which take the form 

C B =^(ny9 M -m)tf, 

~\ E (i^-m^WAW"), (32) 

Ant = " - %, 7/ ^>* 

- fiff^Tii^* - ^e 7 4(l + r 3 )A^, 
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with the field strength tensor Fj$ = d^Au — d v Ay , the meson fields A^\ and the 
electromagnetic field The bold symbols are isovectors, the 7^ are the Dirac 7 matrices 
and \& is a nucleon field which consists of Dirac 4-spinors with isospin space. The masses are 
the baryon mass m = 938.9 MeV and the meson masses m a = 520 Mev, m u = 783 MeV and 
m n = 770 MeV according to a parameter set for the linear model from Horowitz and Serot 
cited as L-HS in [15j. The coupling constants of this parameter set are g a = 10.4814, 



g w = 13.8144 and g p = 8.08488. The charge of the electron e = a/ at He/ Air where a is the 
fine structure constant and he = 197.32 MeV fm. 

Applying the equations of motion and taking the static limit we obtain in the Hartree 
approximation the static Dirac equation ([171]) 

e a ^) a = (ap + V + P(m- S))ip a . (33) 



where a and (3 are matrices like in [18(, e a the single-particle energy of the state ip a , p the 
momentum operator and S and V the the scalar and vector fields 

S = -g^a 

V = gj^ + \g p r z A^ + e |(1 - r 3 )4 7) . (34) 

For the mesons fields we get Klein-Gordan-equations. After neglecting retardation effects 
and taking the Hartree-approximation the meson field equations read 

(-A + ml)<S> a = -g aP s 
(-A + ml)A^ = g uP 
(-A + mJ)^ = \g pP , 

-A A™ = e Pc (35) 

with the scalar density p s , the baryon density p, the The densities are calculated taking 
into account only the occupied positive energy states in the Fermi-sea and neglecting the 
negative energy states in the Dirac-sea ("No-sea" approximation) 

N 



a=l 

N 

a=l 

N 

P3 = ^ Va ^gT 3 J ljj a 
a=l 

N - 1 

^rj a ^ a -{l - r 3 )7o^ Q (-p e ), 



Pc 



a=l 



where r] a are the occupation numbers determined by the BCS-formalism. The electron 
density p e has been considered for the Wigner-Seitz cell calculations but not for studies of 
finite nuclei. 
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B. Solving the Triaxial Dirac Equation 

The solution of the Dirac equation for nucleons in a cubic box differs of course in some 
aspects from that in the spherical one. Therefore we briefly outline the numerical solution 
in the following. 

We decompose ip a in an upper and lower component Pauli spinor: 

*- = fe) <36) 

Hence after applying the transformation e a — > s a — m to energy levels without rest-mass the 
Dirac equation becomes 

e a (p a = erp Xa + U v ip a 

e a Xa = crp tp a + U x Xa (37) 



with the potentials 



u v = -s + v 



U x = -2m + S + V. (38) 

Now we obtain an "effective Schroedinger equation" by inserting the lower component into 
the equation for the upper one. This method is according to Reinhard 15[ the most efficient 
way to solve the Dirac equation. First we modify the lower component 

(e a - U x ) Xa = <?P Pa (39) 
then we introduce an "effective mass term" which depends on the wave function 

Ba = -^yr ( 4 °) 

and finally get the "effective Schroedinger equation" 

£a¥a = crpB a crp ip a + U v ip a . (41) 

So far the procedure corresponds to the method employed in calculations assuming spherical 
symmetry [19|]. Using the discretization in a cartesian box, however, requires a different 
treatment of angular momentum and spin-orbit terms. With the help of the following 
formula for vector fields A and B commuting with cr 

crA crB = AB + icr (A x B) (42) 

we obtain from the relativistic kinetic energy term 

apB a a-pip a = -VB a Vip a - i (VB«) • (V x a) <p a (43) 

which is like the non-relativistic kinetic energy term plus the spin-orbit term for the upper 
component. From a further modification we obtain an expression called "effective Hamilto- 
nian" ready for implementation 

H v , a ^Pa = -B a A (p a (VB a ) ■ (y<p a ) - % (VB a ) ■ ( V x a) <p a + U v <p a . (44) 
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In order to calculate the eigenvalue e a we can't use the "effective Hamiltonian" like a normal 
Hamiltonian because we have to take into account the effects of the lower component. This 
we do in the following manner: 

Xa = B a (Tp Lf a (45) 

and hence the next approximation of the eigenvalue in the iteration scheme gets 

4 n+1) = / d 3 r & a H^ a + e^xlXa) (46) 

This means that in the lower component the whole new information is contained in the new 
Pauli spinor. The total binding energy is calculated like in 20| using cartesian coordinates 



E = y~] 7] a e a - ^ / tVr 

a 



g a $ a {r)p a {r) 



(47) 



+ g u A^p(r) + \g p Af Pz {r) + eA^ Pc (r 
_l_ e + E ■ 

~ -'-'cm ~ -^pair 

with a center of mass correction in the case of finite nuclei 

E cm = -\hjjj with /^ = 41 A~ 1/3 MeV (48) 

in compliance with 21] and a pairing energy E pair described in the next section. 

For the variation of the wave functions in the cubic box we employ once more the imagi- 
nary time step in the following manner: First we operate on the upper component with the 
imaginary time step and the "effective Hamiltonian": 

^ = exp(-A^, Q ) <pM (49) 

then the lower component is calculated: 

X { : +1) =B a ap^ (50) 

and finally both components are orthonormalized together via the Gram-Schmidt method 
considering the symmetries of the Dirac spinors. The symmetries are the same as in 
the Skyrme-Hartree-Fock calculations which are time reversal invariance, parity and z- 
signature. These symmetries furthermore prevent the solution from "slipping" into the 
Dirac sea. In case of a Dirac spinor these symmetries result in parity properties summarised 
in Table [Til which corresponds to the Dirac spinor ansatz in spherical symmetry written in 

3- 





x=0 y=0 z=0 


Re (p a (r, +i) 
Im cp a (r, +\) 
Re ip a {r,-\) 
Im ip a (r, -i) 


+ + Pa 
~ ~ Pa 
- + ~Pa 
+ - ~Pa 


Re Xa(r, +5) 
Im Xa{r,+^) 
Re Xa(r, -5) 
Im Xa(r, -5) 


~ ~ ~Pa 
+ + ~Pa 
+ - Pa 
~ + Pa 



TABLE II: Parity properties of the Dirac spinor with respect to the coordinate planes 
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For a comparison the energy of homogeneous asymmetric nuclear matter is calculated 
similar to [21 



C. Numerical Procedure 

The numerical method for solving the equations for the baryonic wave functions (1491) and 
( 1501) is essentially the same as in the case of the Skyrme Hartree-Fock approach. Thus we 
restrict the discussion in this section to the solution of the meson field equations ( 135|) and 
add some comments on the imaginary time step. 

The meson equations have been solved with a finite difference scheme employing the 
conjugate gradient iterator operating on the meson fields with periodic boundary conditions. 
The conjugate gradient iterator has the numerical advantage that there is no operator matrix 
needed but only the operation of the differential operator on the meson field. The conjugate 



gradient method has been developed to solve linear equations [24], |25j and is now applied 
to a whole variety of numerical problems for example to finite element solver for elliptic 
boundary value problems on an adaptive mesh with hierarchical basis preconditioning [26], 
which provides a very fast algorithm. The main idea of the conjugate gradient step is to 
solve the linear equation Ax — b = with the linear operator A and a vector b by searching 
the minimum of the quadratic form 



q{x) 



\x T Ax-b T x. (51) 



In order to search the solution numerically one can use an iteration scheme following the 
gradient method. Then it was discovered that the iteration is accelerated if one searches not 
straight in gradient direction but in the hyper-plane perpendicular to all previous directions. 
Theoretically the conjugate gradient step converges in less or equal steps than the dimension 
of the vector space. In practical applications the machine errors require a restart after a 
certain amount of steps. 

The overall numerical procedure has a good convergence. In the imaginary time step the 
step At for A = At/H could be set to 4.0 x 10 _24 s which is even larger compared to the 
corresponding Skyrme calculations. In the test runs we obtained results with the parameter 
set L-HS which agree with [16| within numerical accuracy. We used this parameter set also 
for the actual calculations to compare the main properties of the relativistic mean field with 
the Skyrme calculations in the Wigner-Seitz cell. 

IV. PAIRING CORRELATIONS 

Various properties of a neutron star, like e.g. its fluidity, the opacity with respect to 
neutrino propagation etc., are very sensitive to occurrence of pairing correlations. Therefore 
we included the possible effects of pairing in all calculations. Our special attention was 
focussed on isospin T = 1 pairing for nucleon pairs with total momentum equal to zero in 
the 1 Sq partial wave like in an earlier approach in a spherical box [g|. Using the standard 
BCS approach the pairing gap A*, for pair of nucleons with momenta k and — k is obtained 



by solving the gap equation [27 



9 r°° a 

A fc = --/ dk'k' 2 V(k,k') -. k ' . (52) 

Jo 1 ^Jie'-e^ + Ai, 1 J 



7T 
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Here V(k, k') denotes the matrix elements of the NN interaction in the 1 So partial wave, e k 
the single particle energy for a nucleon with momentum k and Bp the Fermi energy. 

Instead of using the matrix elements of a realistic NN-interaction which is fit to the 
scattering data we have decided to use the density dependent zero range effective interaction 



by Bertsch and Esbensen [28 



V{v u r 2 ) =V (l-K f^^j ) 6( ri - r 2 ) (53) 

with the parameters Vq = 481 MeV fm 3 , k = 0.7, a = 0.45 and the cut-off parameter for the 
gap-equation e c = 60 MeV. These parameters were derived from a realistic NN interactions 
by Garrido et al. (29| . 

The occupation probabilities r\ k = v\ which are used to define the densities of the Skyrme 
Hartree-Fock or the relativistic mean field approach are determined from the quasi-particle 
energies E k [9J 

- ifl-fi^t) (54) 



2 V E. 



with 



E k = J(e k -e F y + Al (56) 



in which the pairing gap A k , the single particle energy e k and the Fermi energy Bp enters. 
The BCS-equations have to be solved in a self-consistent procedure fixing the Fermi energy 
ep by the particle number condition for iV nucleons: 



n = J2 v I ( 57 ) 



From the coefficients u k and v k of the standard BCS approach [9j and the corresponding 
single-particle wave functions ip k one can calculate the anomalous density 

X(r) = u k v k \(p k (r)\ 2 . (58) 

k 

For a zero range pairing interaction as the one of eq. (1531) . a a local gap function can be 
defined: 

A(r) = -l/(r) X (r). (59) 
The pairing correlations for continuous asymmetric nuclear matter have been evaluated 



using the techniques described in [27|, [30 



V. RESULTS AND DISCUSSIONS 

In the first part of this section we are going to discuss the results of Hartree-Fock calcula- 
tions using the Skyrme force with the parameter set SLy4 as defined in |l0j . The calculations 
are performed in a Wigner Seitz (WS) cell with a shape of a cubic box. The size of the box 
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x,y,z [fm] x,y,z [fm] 



FIG. 1: (Color online) Density distributions resulting from Skyrme HF calculations for 
protons (black color) and neutrons (red color) as a function of cartesian coordinates x, y, z. 
The panels in the left column refer to densities 0.0166 fm~ 3 (top), 0.0317 fm -3 , and 0.0565 
fm -3 (bottom), while those in the right column are obtained for baryon densities 0.0681 
fm -3 (top), 0.079 fm -3 , and 0.1 fm -3 . Further discussion in the text. 

R has been assumed to be identical in all 3 cartesian directions and has been adjusted to 
minimise the total energy per nucleon for the density under consideration. The calculations 
have been performed for charge neutral matter containing protons, electrons and neutrons 
in /^-equilibrium. 

A few typical density distributions resulting from these variational calculations are dis- 
played in Fig. [T] with densities increasing from top to bottom and all densities displayed in 
the left column being larger than those in the right part of the figure. 

We start our discussion with the top panel in the left column representing a nuclear 
structure at a baryonic density of 0.0166 fm~ 3 . In this case the density profiles are identical 
in all 3 cartesian directions. This means that we obtain a quasi-nuclear structure with 
spherical symmetry in the center of the WS cell. The proton density drops to zero at a 
radial distance of around 4 fm. The neutron density profile drops around the same radius 
from a central density of around 0.1 fm -3 to the peripheral value of around 0.01 fm -3 . This 
means that at this density we have obtained a structure of quasi-nuclear droplets forming a 
cubic lattice, which is embedded in a sea of neutrons. 

The second panel in the left part of Fig. [T] displays the density distributions, which have 
been obtained at a density of 0.0317 fm~ 3 . In this case we obtain deformed quasi-nuclear 
droplets with radii, which are slightly larger in one direction (chosen to be the ^-direction, 
dashed curves) than in the other two, which means that we find prolate deformation. 

At slightly larger densities the deformation of the quasi-nuclear structures increase until 
we reach a density at which the proton density does not vanish along one of the three axis. 
Such an example (baryon density 0.0565 fm -3 ) is displayed in the bottom panel of the left 
column. In this case we have quasi-nuclear structures in the shape of rods parallel to the 
z-axis. The density of these rods is not homogeneous along the symmetry axes. Note that 
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FIG. 3: Profiles for the proton density distribution forming a slab-structure at a density of 

0.0775 fm- 3 . 

in this example size of the WS cell became so small (R=10 fm) that the distance from the 
center to the boundary of the WS box lies within the range displayed in this figure and 
therefore has the boundary been indicated by the dotted line in this panel. This structure 
is also displayed in Fig. [2J where the profile of the proton-density is displayed in the xy and 
xz plane, respectively. 

Performing HF calculations at a density of 0.0681 fm -3 led to a density density distribu- 
tion as displayed in the top panel of the right column in Fig. [TJ In this example the proton 
as well as the neutron density is essentially constant in the (x, y, z = 0) plane. As a function 
of the third coordinate (z, dashed lines) the proton density is reduced from the central value 
at z = to zero at the border of the WS cell and also the neutron density is reduced by 
about 25 percent going from the central to the peripheral values of z. Therefore in this 
case we observe a structure in form of parallel slabs. This slab structure is also displayed in 
Fig. [31 From this presentation in particular it gets obvious that the density within such a 
slab at z = is not really a constant but drops in particular along the diagonals of the WS 
cell with x = y, z = 0. 

At even larger densities the Skyrme Hartree-Fock calculations in a cubic WS cell yield 
structures, with smaller neutron densities in the center of the WS cell as compared to the 
boundaries. An example of such an inverse structure, which corresponds to bubbles in the 
sea of nuclear matter, is displayed in the second panel of the right column of Fig. [1] at a 
density of 0.079 fm -3 . The proton density, which is hardly visible in this example, drops 
from a peripheral value of around 0.004 fm~ 3 to a central value of zero. 

As a final example we present in the bottom panel of the right column of Fig. [T]the results 
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Skyrme 


RMF 




HF 


TF 


H 


TF 


droplet-rod 


0.042 


0.066 


0.070 


0.062 


rod-slab 


0.070 


0.078 






slab-homogeneous 


0.080 


0.085 


0.075 


0.072 



TABLE III: Comparison of densities at which shape transitions occur using the Skyrme 

and Relativistic Mean Field (RMF) approach. Results are compared, employing the 
microscopic Hartree-Fock (HF), Hartree (H) or the Thomas- Fermi (TF) approach. All 

entries are presented in fm -3 . 

of the HF calculation at a baryonic density of 0.1 fm -3 . At this and larger densities, the 
variational calculation yields homogeneous nuclear matter in /3-equilibrium. This example 
also demonstrates that the cartesian box allows for a clean representation of the limit of 
homogeneous matter. This is in contrast to calculations employing a spherical WS cell. 
Depending on the boundary conditions used, calculations within such a spherical box can 
lead to density profiles, which either show a maximum or a minimum at the boundary. Even 
if one tries to use a set of boundary conditions, which minimise this effect, the resulting 
density profile does not correspond to the homogeneous solution (§]. 

From this discussion we see that the HF calculations in a cartesian WS cell for densities 
in the range of 0.01 fm -3 to 0.1 fm -3 leads to quite a variety of shapes and quasi-nuclear 
structures with smooth transitions in between. Following the discussions above these struc- 
tures may be characterised as quasi-nuclei, rod-structures, slab structures (all embedded in 
a sea of neutrons) and, finally, the homogeneous matter. The densities at which the transi- 
tions from one shape to other occur according to our Skyrme HF calculations are listed in 
table HTT1 The transition densities are very similar to those obtained in [3l| . 

The energies per nucleon and the proton abundances resulting from Skyrme Hartree- 
Fock calculations are are displayed in the lower and upper panel of Fig. HI respectively. The 
solid lines indicate the results for the evaluation of homogeneous matter in /3-equilibrium. 
The results of calculations performed in cubic WS cells are presented in terms of individual 
symbols. Those symbols, which scatter around the homogeneous matter results are obtained 
from WS calculations, constraining the HF single-particle wave functions to plane waves. 
Therefore the scattering of these homogeneous matter calculations within WS cells of finite 
size around the homogeneous result for infinite matter is a measure of the shell-effects in 
the WS calculations on the calculated energy and proton abundances. 

The Hartree-Fock calculations, which allow for the formation of inhomogeneous quasi- 
nuclear structures, lead to a reduction of the calculated energy of 1 to 2 MeV per nucleon. 
This gain in energy is reduced with increasing density up to the density of 0.085 fm~ 3 at which 
the energies of the inhomogeneous structures merge into the results for the homogeneous 
matter. At densities below this value of 0.085 fm -3 the balance between the gain in binding 
energy due to a local increase of the baryon density and the loss of binding energy due to 
the localisation of nucleons and surface effects favours the occurrence of inhomogeneities in 
the baryon densities. 

This balance between bulk energy arising from the energy density of nuclear matter 
treated in a local-density approximation and surface effects is also contained in the Thomas- 
Fermi (TF) approach. In this section we want to investigate to which extent the results 
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FIG. 4: (Color online) Proton-abundances and energy per nucleon as obtained from 
Skyrme Hartree-Fock calculations at different densities. The results evaluated in cubic 
Wigner Seitz cells (various symbols) are compared to those of homogeneous infinite matter 
(solid lines) and of Thomas-Fermi calculations. Further details are given in the text. 

of our Hartree-Fock calculations can be reproduced by corresponding TF calculations. For 
that purpose we consider simple parametrisations for the density distribution for protons 
and neutrons, which contain a constant peripheral density p° q ut (q = p or n for protons and 
neutrons, respectively) and an inner part describing the density distribution in the center 
of the WS cell. For spherical quasi-nuclear structures we employ the parametrisation of [H 



Pq{t 



(p. 



pT) 



1 - 



n out 
. rq i 



+ f%*, r<Ri 
R q <r 



(60) 



As an alternative we also consider a Wood-Saxon density parametrisation of the form 



p g (r) = (PT ~ PT) 



1 + exp 



r — r. 



pT 



(61) 



For the description of rod-shape quasi-nuclear structures we use cylindrical coordinates and 
parametrise the dependence of the densities on the radial coordinate in a way correspond- 
ing to eqs. f)60l) or (ISTj) . In the case of quasi- nuclear structures in form of slab-shapes these 
parametrisations are considered for the dependence of the densities on the cartesian coordi- 
nate z. 
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r [fm] 



FIG. 5: (Color online) Density distributions resulting from Skyrme HF calculations for 
protons (black color) and neutrons (blue color) as a function of the distance from the 
center of the Wigner Seitz cell. The densities resulting from HF are compared to those 
determined in Thomas- Fermi (TF) calculations, assuming the parametrisation of f[6"Ul) . 
dotted line, and (IBTl) . dashed line. The example refers to a global baryon density of 0.0166 

fm -3 . 

Assuming those density distributions, the TF energy is calculated as a sum of the bulk- 
energy, i.e. the integrated nuclear-matter energy densities, plus the contribution of a surface 
term of the form [3|, |j] 

E SUTi = F I d 3 r |Vp| 2 . (62) 

^WS-cell 

The parameters of the density distributions in (1601) and (16 ip are varied to minimise the energy 
of the system under consideration. The Parameter F Q for the surface energy term in (1621) has 
been adjusted in two different ways. In a first approach we have considered the properties 
of the nucleus 208 Pb and adjusted F in such a way that the TF calculation reproduced the 
energy and radius of this nucleus derived from Skyrme HF. This leads to a value of Fq of 
68.3 MeV fm 5 and 59.7 MeV fm 5 using the parametrisation of eq. fl60|) and the Wood-Saxon 
parametrisation of eq. flSTl) . respectively. 

Adjusting the surface parameter Fq in this way, one can evaluate the energies of quasi- 
nuclear structures in a WS cell using the TF approximation. The results for these TF 
energies are presented by the dashed dotted line in the lower panel of Fig. HI One finds that 
this procedure leads to energies, which are consistently larger than those obtained in the 
HF calculations. It seems that the TF approach, as it is used here, is underestimating the 
gain in energy due to the formation of inhomogeneous structures. This could be a general 
problem of the TF approximation or a result of the limitation in the variational ansatz for 
the density functions. 

To investigate these possibilities we have considered the different parametrisations dis- 
played in eqs. (!60|) and (ISTj) . It turns out that these two parametrisations lead indeed to 
different density distributions, as displayed in the example of Fig. [5j but it turns out that 
the resulting energy predictions do not exhibit significant differences, so that we present 
only one example for the TF approach in Fig. HI 

We then readjusted the the surface term in (I6"2"|) to obtain an optimal fit of the HF energies 
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for the quasi- nuclear structures in /^-equilibrium. This readjustment of the surface term leads 
to values of the surface parameter F , which are about a factor of one half smaller than 
obtained from the fit to the properties of 208 Ph. Using these readjusted surface parameter 
we observe critical densities for the shape transitions of the quasi-nuclear structures from 
droplets to rods to slabs and to homogeneous nuclear matter at values which are similar to 
the results obtained in the HF calculations. If, however, one uses this reduced values derived 
from the fit to inhomogeneous matter in /3-equilibrium the TF calculation do not give an 
accurate description of Hartree-Fock energies, in which the proton abundance has been fixed 
e.g. to a value of 10 percent. This result can be taken as an indication that in addition to 
the iso-scalar surface term of (162]) an iso-vector surface term might be required in addition 
to obtain a reliable TF approximation to the results of corresponding HF calculations over 
a wide range of proton-neutron asymmetries. 

The upper panel of Fig. H] contains results on the proton abundances for baryonic matter 
plus electrons in /3-equilibrium. The value of the proton abundance assuming homogeneous 
matter increases with density reaching a value of about 4 percent at a baryonic density of 0.1 
fm -3 . Allowing for inhomogeneous, however, this value is almost constant around 3.2 percent 
in the density interval from 0.03 to 0.08 fm~ 3 and yields even larger values for densities below 
0.03fm~ 3 . This trend is also reproduced in the TF calculations. The increase of the proton 
abundances at smaller global densities reflects the fact that at those small densities we 
observe local structures in the center of the WS cells, with large local densities. The proton 
abundance in these quasi-nuclear droplets is significantly larger than the proton abundance 
in the homogeneous matter with the same global density. The scattering of the results for 
the proton abundances as a function of density resulting from the HF calculations reflects 
the shell-effects, which preferentially yield quasi-nuclear with closed shells for the protons. 

A comparison of energies resulting from relativistic mean field calculations in a Wigner 
Seitz cell are displayed in the lower panel of Fig. [61 Comparing these results with the 
corresponding values displayed in Fig. H] one finds that the energy gain due to the formation 
of inhomogeneous structures is much weaker in the relativistic mean field calculations as 
compared to the Skyrme model. This is also reflected in the corresponding Thomas-Fermi 
calculations. Note that also in this case we have adjusted the constant F of the surface 
term in (I6"2"l) to reproduce the bulk properties of 208 Pb as predicted by the relativistic mean 
field calculations. This leads to value for F of 87.4 MeV fm 5 and 80.3 MeV fm 5 using 
the parametrisation of eq. (l6"Ul) and the Wood-Saxon parametrisation of eq. flBT]) . respectively. 
Both values are significantly larger than the values required for Fq in the case of the Skyrme 
model used above. 

The different interplay between volume-, surface-, symmetry- and Coulomb effects in the 
relativistic mean field model as compared to the Skyrme model also leads to smaller values for 
the proton abundance in the region of nuclear densities, in which inhomogeneous structures 
emerge. The values around p = 0.02 fm— 3, displayed in the upper panel of Fig. El are about 
40 percent smaller than the corresponding values obtained in the Skyrme model (see Fig. Hj). 
The differences in the balance between volume- and surface-contributions to the energy also 
lead to different quasi-nuclear structures in the nuclear models under consideration. It is 
worth mentioning that within the relativistic mean field mode we do not find any formation 
of slab-like structures. Therefore the table HTT1 contains for this case only transition densities 
for droplet to rod structures and the formation of a homogeneous structure. 

The density profiles obtained from these 2 approaches also yield different results. As an 
example we present in Fig. [7] the density profiles at p = 0.032 fm -3 , a density at which both 
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FIG. 6: (Color online) Proton-abundances and energy per nucleon as obtained from 
relativistic mean-field calculations at different densities. The results evaluated in cubic 
Wigner Seitz cells (various symbols) are compared to those of homogeneous infinite matter 
(solid lines) and of Thomas-Fermi calculations. Further details are given in the text. 




FIG. 7: (Color online) Density profiles for protons and neutrons as derived from Skyrme 
HF and relativistic mean field calculations at a global density of p = 0.032 fm -3 . 




FIG. 8: (Color online) Local pairing gap A(r) (upper panel, see eq. (l59|) ) and anomalous 
density x( r ) (lower panel, see eq. (l5Bl ) for the configurations, which are also considered in 

Fig. m 

the relativistic as well as the Skyrme model yield a droplet structure. Note, that in the case 
of the Skyrme calculation we obtain a Wigner Seitz cell with a length of 26.4 fm which leads 
to a borderline as indicated by the dotted line, while the corresponding borderline for the 
RMF calculation is identical to the frame of the figure. 

Finally, a feature of the pairing correlations obtained in these calculations shall be dis- 
cussed. For that purpose we present in the upper panel of Fig. [8] the local pairing gap A(r) 
(see eq. (1591) for the formation of neutron pairs, as obtained in the Skyrme and relativistic 
mean field model at p = 0.032 fm -3 . In both of these approaches one observes a suppression 
of the local gap A(r) in the region of the quasi-nuclear structure, i.e. in the region where 
the density is large. 

This phenomenon has been observed before [8, 31, 32|, 33] and has lead to discussions about 



various phenomena, which are related to to this periodic structure of the gap parameter. It 
should be noted, however, that this suppression of the gap parameter in the high-density 
region of the quasi-nuclear structure is either to the local-density approximation, which 
is used to calculate this local gap or to the assumption of the density-dependence of the 
interaction strength for the pairing interaction, like the one, which we have considered in 
our calculations (see eq. (l53l ). If, rather than looking at the local gap parameter A(r), we 
inspect the anomalous density x{ r ) ( see eq. (l58l) ). one finds even a small enhancement of 
the anomalous density in the region of the quasi-nuclear structure. This suggests that the 
reduction of the pairing gap in the region of high densities might be an artefact of the special 
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interaction considered. 



VI. CONCLUSIONS 

The structure of neutral baryonic matter is investigated in a region of baryon densities 
between 0.01 and 0.1 fm -3 performing various Hartree-Fock an mean-field calculations with 
inclusion of pairing correlations in a periodic lattice of Wigner-Seitz (WS) cells of cubic 
shapes. In this region of densities, which should occur in the crust of neutron stars, one 
observes structures ranging from neutron-rich nuclei embedded in a sea of neutrons up 
to homogeneous matter. The symmetries of the WS cell allow the formation of triaxial 
structures but also include rod- and slab-like structures and provide a natural transition to 
the description of homogeneous matter. 

For the baryonic components a Skyrme Hartree-Fock approximation has been considered 
as well as a relativistic mean field model. Both approaches yield an intriguing variety 
of quasi-nuclear structures with smooth transitions in between. The occurrence of special 
structures as well as the critical densities at which transitions between those structures occur 
depend on the nuclear model considered. 

The resulting energies as well as the proton abundances can fairly well be reproduced by 
a Thomas- Fermi approach, if the constant, determining the strength of the surface term is 
adjusted to reproduce the results of the microscopic calculations. A surface term depending 
on the isospin asymmetry might be required to obtain Thomas-Fermi results, which are 
reliable over a large interval of proton-neutron asymmetries. 

Pairing-correlations have been evaluated within the BCS approach, assuming a density- 
dependent contact interaction. This leads to local pairing gaps for neutron pairing, which 
are significantly smaller in the regions of the quasi-nuclear structures as compared to the 
bulk of the neutron sea. It is argued, however, that this feature might be an artefact of the 
density-dependence of the effective pairing interaction. 

The present studies provide an interesting starting point for further studies on the prop- 
erties of matter in the crust of neutron stars. The single-particle energies and wave-functions 
could be used for a microscopic study of response-functions, which allow e.g. the evaluation 
of neutrino opacities. 

This has been supported by the European Graduate School "Hadrons in Vacuum in 
Nuclei and Stars" (Basel, Graz, Tubingen), which obtains financial support by the DFG. 



C.J. Pethick and D.G. Ravenhall, Ann. Rec. Nucl. Part. Sci. 45, 429 (1995). 
M. Hashimoto, H. Seki, and M. Yamada, Prog. Theor. Phys, 71, 320 (1984). 
K. Oyamatsu, Nucl. Phys. A 561, 431 (1993). 
K. Oyamatsu and K. Iida, nucl-th/0609040 



T.H.R. Skyrme, Nucl. Phys. 9, 615 (1959). 

D. Vautherin and D.M. Brink, Phys. Rev. C5, 626, (1972). 
P. Bonche and D. Vautherin, Nucl Phys. A 372, 496 (1981). 

F. Montani, C. May, and H. Miither, Phys. Rev. C69, 065801 (2004). 

P. Ring and P. Schuck, The Nuclear Many Body Problem, (Springer, New York, 1980). 

E. Chabanat, P. Bonche, P. Haensel, J. Meyer, R. Schaeffer, Nucl. Phys. A635, 231, (1998). 



23 



P. Bonche, H. Flocard, P.-H. Heenen, S.J. Krieger, M.S. Weiss, Nucl. Phys. A443, 39 (1985). 
K.T.R. Davies, H. Flocard, S. Krieger, and M.S. Weiss, Nucl. Phys. A342, 111 (1980). 

D. Vautherin, Phys. Rev. C7, 296 (1973). 

J.D.Jackson, Classical Electrodynamics (2nd ed., Wiley, New York 1975). 

P.-G. Reinhard, Rep. Prog. Phys. 52, 439 (1989). 

C.J. Horowitz and B.D. Serot, Nucl. Phys. A368, 503 (1981). 

R. Fritz, Korrelationen una 1 Relativistische Effekte in Atomkernen (Ph-D thesis, Eberhard- 
Karls-Universitat Tubingen, 1994). 

J.D. Bjorken, S.D. Drell, Relativistic Quantum Mechanics (McGraw Hill, New York, 1964). 
M. Rufa, P.-G. Reinhardt, J. Maruhn, W. Greiner, and M.R. Strayer, Phys. Rev. C38, 390 
(1988). 

B.D. Serot and J.D. Walecka, Adv. Nucl. Phys. 16 (1986) 

F. Hofmann, CM. Keil, and H. Lenske, Phys. Rev. C64, 034314 (2001) 

M. Kleinmann, R. Fritz, H. Miither, and A. Ramos, Nucl. Phys. A579, 85 (1994). 

S.S. Avancini, M.E. Bracco, M. Chiapparini, and D.P. Menezes, Phys. Rev. C67, 024301, 

(2003). 

M.R. Hestenes, and E. Stiefel, Nat. Bur. Standards, J. of Res. 49, 409 (1952). 

J.K. Reid, Large Sparse Sets of Linear Equations, 231 (Academic Press, London and New 

York, 1971). 

H. Yserentant, Appl. Math. Comput. 19 no. 1-4, 347 (1986) 

J. Kuckei, F. Montani, H. Miither, and A. Sedrakian, Nucl. Phys. A723, 32 (2003). 

G. F. Bertsch and H. Esbensen, Ann. Phys. 209, 327 (1991). 

E. Garrido, P. Sarriguren, E. Moya de Guerra, and P. Schuck, Phys. Rev. C60, 064312 (1999). 
S.A.Fayans, S.V. Tolokonnikov, E.L. Trykov, and D. Zawischa, Nucl. Phys. A676, 49 (2000). 
P. Magierski and P.-H. Heenen, Phys. Rev. C65, 045804 (2002). 

P. Magierski, Phys. Rev. C75, 012803 (2007). 



M. Baldo, E.E. Saperstein, and S.V. Tolokonnikov, nucl-th/0609031 



